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SUMMARY 

This paper presents the application of a class of multi-grid methods to 
the solution of the Navi er-Stokes equations for two-dimensional laminar flow 
problems. The methods consist of combining the full approximation scheme - 
full multi-grid technique ( FAS-FMG) with point-, line-, or plane-relaxation 
routines for solving the Navi er-Stokes equations in primitive variables. The 
performance of the multi-grid methods is compared to that of several single- 
grid methods. The results show that much faster convergence can be procured 
through the use of the multi-grid approach than through the various sugges- 
tions for improving single-grid methods. The importance of the choice of 
relaxation scheme for the multi-grid method is illustrated. 


INTRODUCTION 

In spite of the rapid development of computer hardware over the last two 
decades, there is need for continuous improvement of computational techniques 
in order to achieve better accuracy with minimal computational effort. In 
most numerical methods, for fluid flow calculations, better accuracy is 
obtained through the use of higher-order discretizations, or finer grid distri- 
butions, or a combination of both. With a large number of grid points, direct 
methods of solution are often not feasible, because of memory limitations, so 
an iterative method must be used. However, it has been found that Increasing 
the number of grid points distributed in the computational domain leads to a 
deterioration of the convergence rate. Thus, for most practical problems in 
two- or three-dimensions, the grid distribution required for accurate solution 
would lead to excessive computational times (of the order of hours for two- 
dimensional and days for three-dimensional calculations), even on very large 
mainframe computers. Hence, the quest for more efficient solution procedures. 
The multi-grid technique is emerging as a very promising tool for accelerating 
the convergence rate of iterative procedures, especially for calculations with 
very fine distributions. 


*This paper is reproduced from the Proceedings of the Institution of 
Mechanical Engineers, Part C, 1989, by permission of the Council of the 
Institution. 

*Work funded under Space Act Agreement C99066G. 



The technique for accelerating convergence of an Iterative procedure 
through the use of multiple grids was applied as early as 1935 by Southwell 
(ref. 1). Other works along similar lines were reported by Southwell (refs. 2 
and 3), Stiefel (ref. 4), Federenko (ref. 5) and Wachspress (ref. 6), among 
others. These were all two-grid level methods. The idea of utilizing fully 
multiple-grid levels was introduced by Federenko (ref. 7) for a Poisson-type 
problem on a rectangular gird, and the approach was generalized by Bakhalov 
(ref. 8) to any second-order elliptic operator with continuous coefficients. 
The first actual numerical computations with the full multi-grid technique 
were reported by Brandt (ref. 9) for some boundary-value problems, using a 
finite-difference approach. The multi-grid method has been applied with 
finite-element formulations by N 1 col aides (ref. 10). In most of these, linear 
differential equations were solved. More recently, multi-grid techniques for 
solving nonlinear equations have been reported by Brandt (ref. 11), Ghia 
et al . (ref. 12) and Vanka (ref. 13). 

The previous studies show that the multi-grid methods are quite promising 
but a lot of work still needs to be carried out in order to realize their full 
potential. In the present study, the multi-grid technique is applied, with 
some finite-difference methods, to solve the Navier-Stokes equations for two- 
dimensional model flow problems. Their performance is compared with that of 
various single-grid methods for varying grid distributions and Reynolds 
numbers . 


MATHEMATICAL MODEL 
The Multi-grid Concept 

Consider a differential problem represented over a given domain D by the 
linear equation: 


L U(x) = F ( x ) (1 ) 

where L is the differential operator, U is the unknown variable, F is a 
known function, and x = (x 1( x 2 , ... x d ) are the d independent variables of 
the d-dimensional problem. If the domain is divided Into a computational 
grid with spacing h, the finite-difference equivalent of the equations may be 
written as: 


Lhyh( x h) _ f* 1 (x* 1 ) (2) 

Solution of this equation by an iterative procedure such as Gauss-Seidel , 
Jacobi, line relaxation, etc., have been found to converge rapidly only for the 
first few iterations, and more slowly thereafter. By considering a Fourier 
analysis of the error-reduction process of typical relaxation procedures, 

Brandt (ref. 11) showed that they are only efficient in smoothing out those 
error components whose wavelengths are comparable to the grid mesh size. 

Error components with longer wavelengths are smoothed out at comparatively 
slower rates. Thus, the idea of the multi-grid technique is to smooth out 
high frequency error components by performing a few iterations of the relaxa- 
tion process on the fine grid. The remaining errors are then transferred to a 
coarser grid where the corresponding lower frequency error components are again 
smoothed out in a few iterations. Further transfers are made to even coarser 
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grids and the process Is repeated until all the error components have been 
smoothed out. The results are then progressively transmitted back to the 
finer grids. The overall effect of this procedure Is that the various Fourier 
components of the error are removed on grid meshes most efficient for the pur- 
pose, thereby accelerating the convergence rate on the fine grid on which the 
solution to the finite difference equation Is sought. 

Thus, in the multi-grid method, the computation Is carried out on a 
series of grid meshes D k with the corresponding grid functions U k , where 
k = 1,2,3, ... M, with k = M representing the finest mesh, and the meshes 
becoming coarser the lower the value of k. The mesh sizes usually differ by 
a factor of 2, so that hk+i = 1/2 h^. Although this is not a requirement of 
the multi-grid method, it enables simple weighting functions to be employed In 
Interpolation routines. The finite difference equation can then be written 
for the k^ grid as: 


L k U k = F k 


(3) 


There are several algorithms for implementing the multi-grid idea, each 
with several possible variations. One of the simplest is the correction scheme 
(CS) which is applicable to linear problems. In this scheme, the calculation 
starts on the finest grid D M , and an approximate solution to equation (3) is 
computed by a relaxation method such as the Gauss-Sei del , line relaxation, etc. 
Unless the approximate solution U k (fortuitously) satisfies the difference 
equation (3) and the boundary conditions, there will be a residual R k given 
by: 


F k - L k U k = R k 


(4) 


A few iterations are performed on grid until the rate of reduction 
of the residuals falls below a desired theoretical rate. The residuals are 
then transferred by restriction to the next coarser grid, D M_ ^ and a correc- 
tion function 6U k_1 is obtained by solving the equation: 

L k-l 5lJ k-l = j k -l R k (5) 


where L k- 1 is the difference operator on grid D k_ l and l£ ' is the 

restriction operator. If k-1 is the coarsest grid, equation (5) is solved 
completely on this grid, otherwise, the problem is transferred to the next 
coarser grid as soon as the rate of smoothing of the residuals falls below the 
required norm. Once equation (5) has been solved, the correction function 
SU k ~l is prolongated to grid D k , and U k is subsequently corrected as: 


U k . U k ♦ £ ,5U k -' 
new old k-1 


( 6 ) 


where l£_^ is the prolongation operator. Similar equations are employed to 

transfer (prolongate) the intermediate corrections from the coarser grids, 
once the prescribed convergence criteria have been satisfied. This process of 
relaxation, restriction and prolongation is repeated until the desired accuracy 
on the finest grid, is achieved. 
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The scheme described above serves to Illustrate the multi-grid Idea, but 
Is Inadequate for nonlinear problems such as the Navler-Stokes equation, of 
interest here. For this, the full approximation storage (FAS) scheme Is 
employed. 

In the FAS scheme, the full approximation Is stored on a coarse 

grid rather than Its correction SU^-I as In the CS scheme. Is given 

by: 


U k 1 = I k 1 U k + SU k_1 (7) 

where Is the fine grid approximation. The full approximation can be 

obtained on the coarse grid by solving the corresponding FAS equation: 

L k-1 (J k-l = L k-l^k-l y k^ + - L k U k ) (8) 

The fine grid approximation is then corrected as: 

U k - U k + I k fu k - ] - I k- ^ U k ^ (9) 

u new ~ U old + A-A U A k u oldy 


Comparison of equations (5) and (8) shows that the CS and FAS schemes are 
equivalent for linear problems. However, equation <8) is additionally suitable 
for nonlinear problems, for which the former Is not valid. 


Present Implementation 

In the present paper, the FAS-FMG (full approximation storage - full 
multi-grid) algorithm originally developed by Brandt (ref. 11) Is employed. 

The procedure commences by solving iteratively until convergence, the finite 
difference equations on the coarsest grid, using a relaxation routine. The 
solution Is then prolongated to the next finer grid, where it serves as the 
Initial approximation of the solution to the difference equations on this grid. 
Improved solutions are obtained by applying the relaxation routine on this 
grid, until the smoothing rate falls below a desired value. The optimum value 
was found by numerical experiments to be a rate of reduction of the normalized 
average residual of 0.8 ( 1 . e . , R n+ i/R n < 0.8). Otherwise, the task of smooth- 
ing the residuals is transferred to a coarser grid and the FAS scheme already 
described is employed. This rate happens to correspond to the theoretical 
smoothing rate computed by Brandt (ref. 11) for the Navier-Stokes equation. 

The iteration on each grid level is continued until the prescribed convergence 
criterion is satisfied, at which time, the solution is transferred to the next 
finer grid level and the whole process repeated. The solution process is ter- 
minated on satisfying this criterion at the finest grid level D^. 

It is only necessary to iterate until full convergence on the finest and 
coarsest grids. Otherwise, the finer grid residuals are considered suffici- 
ently smoothed once the corresponding coarse grid errors have been reduced to 
about 20 percent of the restricted values. Then the corrections are prolon- 
gated back to the fine grid. In the present work, fully adaptive, full multi- 
grid cycles are employed, which are controlled solely by the aforementioned 
convergence criteria and the desired smoothing rate. There is, however, a 
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wide range of possible multi-grid cycles, some of which are discussed by 
Stuben and Trottenberg (ref. 14). 


Restriction and Prolongation 

Restriction and prolongation routines are required to transfer the fine 
grid approximations and residuals onto the coarse grid, and the coarse grid 
corrections onto the fine grid, respectively. The operators were denoted by 

ifc- 1 and l[_ v respectively. In the present work, a staggered grid is 

employed for the velocity components and the pressure, as shown in figure 1. 

So, the transfer operators would be different for each of the variables U,V,P 
and for the residuals Rn and R v . The restrictions are made by averaging the 
neighboring fine grid values, and the prolongations by applying bilinear inter- 
polation to the closest set of coarse grid values. 


Relaxation Methods 

The Navler-Stokes equations can be written for two-dimensional steady, 
laminar plane flow as: 


3(UU) 3(UV) 

1 3P j 

( 3 2 U 

2 ' 
3^U 

3x + 3y " 

p 3x + H 

Ux 2 3y 2 

3(UV) 3(VV) 

1 3P 1 

( 3 2 V 

2 ' 

3^V 

3x + 3y 

> 3y + H 

Ux 2 + 3y 2 , 


( 10 ) 


( 11 ) 


and the continuity equation as: 


3U 9V n 

3x + * 0 


( 12 ) 


Using the staggered grid system (fig. 1), these equations can be discre- 
tized for the 1,j“ cell giving, in each case: 


P - P 

c u 1,j ' 'V'l + I.j T "w u 1-l,j T 'V'l.j+l T "s u 1,j-l + pAx i _ 1/2<j 


a“u, , = aHu 


+ A U U , 


+ A^U, 


+ A U U, 


(13) 


P - P 

+ LLlrJ -LI 

c v i,j - "e’l+I.J T "wM-l.j T "n’l.J+l T "s’i.j-l + P^y 1 j_i /2 


a!v 


= aV 


+ A V V. 


+ A V 


+ aV 


(14) 


!!l±L J ~ U i, j + v < ,3-1 ~ V 1,3 _ o 


(15) 


where the coefficients A c , A e , A n , etc., are derived using the hybrid 
(central /upwi nd) difference scheme proposed by Spalding (ref. 15). The actual 
expressions can be found in Vanka (ref. 13). 
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Similar expressions can be written for the other face velocities by 
Increasing the Indices of 1 and j by one In equations (13) and (14), 
respectively. Thus: 

- P, 


LJ _1±U 

3 


• A e U 1+2,J * A w u U * A n u U.,M * "Ku-l * 

P “ P 

A^V - A^V + A V V + A V V + A V V + — ^ >3 + 1 

A c V 1,j+l “ A e V 1+l , j + 1 + Vl-1,3+1 A n V 1 , j+2 Vl.j pAy 1J+1/2 


( 16 ) 


(17) 


We now assume that the velocities on the right-hand side of equations 
(13), (14), (16), and (17) are known, so that the face velocities (on the left- 
hand side) can be evaluated once the pressures are found. The method for find- 
ing the pressures characterizes the relaxation scheme. In equations (13) to 
(17), there are four unknown velocities and five unknown pressures. Each of 
these can be split into known and unknown components, consisting of guesses 


and 


. Thus, 



U i .3 “ U * 

* u ij 


- 11* 

1 + 1, J " 1+1, j 


> (18) 

V i ,3 = V U 

* v i.j 


- V* 

1 , j+l ‘ V i , j+1 



Vr p u 

+ p ij 

(19a) 

hi, j a p uu 

* P U1.J 

(19b) 

n * 

l-i, j = l-l, j 

* P i-U 

(19c) 

1,3+1 = P i,j+i 

*■ p i.j,l 

( 1 9d) 

- p* 

1 ,j-l “ 1 ,3-1 

* p i .J-i 

(19e) 


We now need to solve the system of equations (13) to (17) for the correc- 
tions, 1 .e. , the ( ' ) ' s 


Plane Relaxation (SIM) 

In equations (13), (14), (16), and (17), each of the velocity corrections 
Is expressed solely in terms of known velocities and pressures and unknown 
pressure corrections, so that we can substitute these equations into equation 
(15) to yield the pressure correction equation. This equation will have five 
unknowns, and it will be necessary to generate the equations at all the inter- 
nal nodes in the two-dimensional plane in order to form a closed set. 
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The penta-dl agonal matrix equation set has the form: 

e, ,p: . = a. ,p' , + b, .p; , , + c, ,p; , i + o, ,p; , , + f, . 

1,3 1,3 1,3 1+1,3 1.3 1-1, j 1,3 1.J+1 1.3 1.3-1 1.3 

where 

A i , j = 1 / ( pA c i + 1 j AX 1,J Ax 1+(l/2),j) 


( 20 ) 


1 .3 


pAc i.j Ax i, 3 Ax i-(i/2),j, 


C l,3 = /V^i.j+l Ayi ’J +(1/2) y 


and 
F 


D 1,j = 7 ^ pA c 1tj Ay 1,j Ay 1.j-(l/2) 

E 1.j = A i,j + B 1,j + C 1.3 + °1 .3 


1.3 


- R U*7 ax l i ft u/ A c 


- R 


1 .3 


1 + 1 


+ 1 


A*,, 




The R's represent the residuals of equations (13) to (17). The equa- 
tion set must be solved simultaneously for all Internal grid nodes in the two- 
dimensional plane. The strongly implicit method (SIM) of Stone (ref. 16) is 
applied with an Iteration parameter of 0.9 for the solution. Once all the 
pressure corrections are obtained, the velocity corrections can be computed 
from: 


u i.J ■ 

K 3 * ( 

P ' 

i-i.j 

- P i.j)/ 

^ ax 1-(l/2>,j ] 1 

1*1 . 
1 ,3 

(21) 

1+1,3 

Km * 

( P U 

N 

- P ' 

*1+1, 3> 

l/ p 4x H<l/2)j] 

K.,., 

(22) 

V i,j = 

Kj * ( 

p 1 

1.3-1 

- p ij )/ 

r p Ay i ,j-d/2>]y 

i ,3 

(23) 

,3+1 = 

Km * 

( P i.3 ■ 

' p i . j+i) 

I'Hl.j.U /2>] 

K.„, 

(24) 
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Due to the nonlinearity of the momentum equations It may be necessary to 
under-relax the velocity corrections so that only a fraction of them Is uti- 
lized In equation (18). The usual practice of Implicit under-relaxation Is 
employed here. The optimum under-relaxation factor was found to be 0.9. 


Line Relaxation (CLSOR) 


Line relaxation Is obtained, using a coupled, line successive over- 
relaxation scheme. If we presume that only the pressure corrections for nodes 
lying along a line (In the x- or y-dlrectlon) are unknown, all others are 
known, presumably zero. Then there will be three unknown pressure corrections 
and four unknown velocity corrections to be determined at each node. Equa- 
tion (20) then reduces, for x-llne relaxation, to one of the form: 


E 1 , j P i , j " A 1,j P i + 1.J + B i,j P i-U + F 1 .3 


(25) 


Similar equations need to be written for the other Internal nodes along the 
line. In order to form a closed set. The system of tri-diagonal equations is 
then solved simultaneously, with a tri-diagonal matrix algorithm (TDMA), which 
is efficient and has a low operation count. The velocity corrections are sub- 
sequently computed, for all the nodes lying along the line, using a variant 

of equations (21) to (24) in which the pressure corrections p \ t y\ and 


p! have been set equal to zero. The procedure Is then repeated for the 

otlier parallel lines In the computational domain. Y-llne relaxation can be 
performed in a similar manner. The optimal under-relaxation factor was found 
to be 0.6 for test case 1 and 0.9 for test case 2. 


Point Relaxation 

This is the so-called symmetric coupled Gauss Seidel (SCGS) method pro- 
posed by Vanka (ref. 13). It assumes that only the pressure correction at the 
central node Is unknown so that there are five equations and five unknowns 
(one pressure correction and four velocity corrections) to be determined at 
each node. These equations can be solved algebraically, since simple manipula- 
tions yield an explicit expression for the pressure correction. The velocity 
corrections are then computed from a variant of equations (21) to (24) in 
which all pressure corrections other than the central one are set to zero, 
the procedure is repeated for all internal grid nodes. A lexicographic visit- 
ing order, which is reversed in alternate iterations, was employed in the 
present study. The optimum under-relaxation factor was found to be between 0.5 
and 0.6 for this method. The SCGS differs from the SIVA (simultaneous variable 
adjustment) method of Caretto et al . (ref. 17) only through the incorporation 
of additional unknown velocities on the right-hand side in the latter. For 
example, on the right-hand side of equations (13) and (14), respectively, 

Lh.i i and V^ , are also considered unknown. Similarly, in equations (16) 
and u7) L), j and are considered unknown. The explicit expression for 

the centralpressure correction is more complex in SIVA. 

Comparison of Relaxation Methods 

In the present relaxation methods, the pressure-velocity coupling, inher- 
ent in the momentum and continuity equations, is retained. Further, all the 
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velocity components are obtained simultaneously. They differ from the SIMPLE 
method of Patankar and Spalding (ref. 18), and proposed derivatives (refs. 19 
and 20) which solve for the velocity components in a sequential or decoupled 
manner. 

The plane-relaxation method is the most implicit, and would normally be 
expected to be the most stable, having also the fastest convergence rate. How- 
ever, the coefficients of the equations require two-dimensional storage arrays, 
so that the computer memory requirements are high. On the other hand, the 
point-relaxation method is the least implicit, exhibiting a space decoupling 
from neighboring grid points. The decoupling effect is reduced somewhat by 
computing, at each node, all the four velocities on the faces of the control 
volume. Hence, in a complete sweep of the computational domain, each velocity 
is updated twice. However, no arrays are required to store any coefficient of 
the equations. The line-relaxation method lies in between. Only the veloci- 
ties normal to the relaxation line are computed twice. Also, the coefficients 
only need be stored along one line of grid nodes. The increase in memory 
requirements is thus minimal. 

The expected superiority of plane relaxation is not fully realized in the 
present study, since the methods are applied to nonlinear problems. The equa- 
tions are linearized and the coefficients have to be repeatedly updated as soon 
as better values are available. The point-relaxation method ensures the most 
rapid updating of coefficients, and the plane-relaxation method, the slowest. 
Hence, the advantages of increased implicitness in the latter are offset, some- 
what, by slower updating of coefficients, and vice versa. 


RESULTS AND DISCUSSION 
Test Cases 

Three test cases were chosen in order to evaluate the performance of the 
present methods. The first is the laminar flow in a lid-driven square cavity 
(see fig. 2(a)). The second case is the laminar flow through a backward facing 
step with an expansion ratio of 2 (see fig. 2(b)). Finally, the lid-driven 
cavity flow was computed in several cavities over a range of aspect ratios. 

For case 1, computations were performed at two Reynolds numbers, namely 
100 and 1000, based on the velocity of the moving lid and the cavity width. 

For case 2, computations were performed at Reynolds numbers (based on the mean 
inlet velocity and channel height) of 100 and 500 and aspect ratios L/H of 
3 and 6, respectively. The longer aspect ratio of the computational domain was 
due to the correspondingly longer recirculation length encountered in the 
latter. 

The final case compares the performance of the multi-grid methods in lid- 
driven cavities with aspect ratios (B/H) ranging from 1/6 to 6, each at a 
Reynolds number of 100. The aspect ratio was varied by changing H, while 
keeping all other parameters unchanged. 

In all these test cases, the same number of grid nodes was utilized is 
both the x- and y-di rections , so that the grid cell aspect ratio would be the 
same as the aspect ratio of the channel or cavity. 


9 



Comparison of Single-Grid and Multi-grid Methods 

In order to evaluate the performance of the present multi-grid methods, 
their convergence characteristics are compared with those of various single- 
grid methods. Three additional single-grid methods are employed, namely, the 
SIMPLE algorithm of Patankar and Spalding (ref. 18), an Improved variant by 
Issa (ref. 19) named PISO, and another by Van Doormaal and Ralthby (ref. 20) 
called SIMPLEC. Full details of these methods can be obtained from the origi- 
nal references and will not be repeated here. It suffices to mention that they 
are all based on the line-by-line solution of the Navler-Stokes equations In a 
decoupled manner. Thus, they differ from the present line-relaxation method 
which solves the Navler-Stokes equations, symmetrically along a line, and In a 
coupled manner. 

Table I shows the comparison of the CPU times and the number of Iterations 
required for convergence when the six single-grid methods and three multi-grid 
methods were applied to the lid-driven square cavity flow. The results are 
presented for six different grid levels. The convergence criterion In each 

case Is that the residual norm R < 10 . R Is defined as: 


R = 


R 5 * R v ^ R m 


3 x (NI-2) x (NJ-2) 


where 

Ru sum of residuals of the U-momentum equation at all modes 

r v sum of residual of the V-momentum equation at all nodes 

R m sum of residual mass sources at all nodes 

(NI-2) number of Internal grid nodes In the x-dlrectlon 

(NO-2) number of Internal grid nodes In the y-dlrectlon 

At this convergence level all methods gave the same results. In the columns 
for the multi-grid methods, three values are presented. The first value repre- 
sents CPU times on the Siemens 7881 Computer. Although some of the computa- 
tions were performed on the Deflnicon DS 1/7 80 processor, test computations 
showed that these take 25 times longer, so the times have been duly converted. 
The second value represents the number of iterations on the finest grid and 
the third gives the total work units which Is the equivalent number of finest 
grid Iterations of the total relaxation work performed on all the levels. The 
latter is thus comparable to the number of iterations in the single-grid 
col umns . 

The results show that among the single-grid methods which solve the equa- 
tions in a coupled manner, plane relaxation produced convergence in the least 
number of iterations, followed by line relaxation, and then point relaxation. 
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However, It required the largest CPU time per iteration, hence, total CPU 
times are comparable for all methods. Contrary to expectations, the decoupled 
single-grid methods generally resulted In faster convergence than the coupled 
ones. This is probably due to the faster rate of updating the coefficients In 
the former, which appears to overweigh the disadvantages of solving the equa- 
tions sequentially. This may not apply to the multi-grid methods, since in the 
decoupled procedure, the residuals, which have to be smoothed on the coarser 
grids, are changed by the pressure correction equation in an inconsistent 
manner . 

Arakawa et al . (ref. 21) reported slight advantages for the coupled, 
multi-grid methods over the decoupled ones. In general, the multi-grid meth- 
ods converged faster than single-grid ones, especially for finer grids. 

Figures 3 and 4 present log-log plots of the CPU times for convergence versus 
the total number of internal grid points. The slopes roughly approach 2 for 
single-grid methods, but are approximately 1 for multi-grid methods, which 
shows that the CPU times for convergence Increases almost quadrati cal ly with 
the number of grid points in the latter, but only linearly in the former. It 
is this independence of the computational work per grid point, on the level of 
grid refinement that makes multi-grid methods quite attractive for fine-grid 
calculations. Usually, less than 30 iterations were required on the finest 
grid. However, in contrast to the single-grid methods, the multi-grid methods 
showed slower convergence with increase in Reynolds number, in the lid-driven 
square cavity flow. Typically, 1.5 to 2.5 times more computational work was 
required at Re = 1000, as compared with Re = 100. This coincides with the 
findings of Vanka (ref. 13) and Barcus et al . (ref. 22). 

At higher Reynolds numbers, the coarser grids appear to become ineffi- 
cient in smoothing out finer grid errors. The point-relaxation method pro- 
duced the fastest convergence of all multi-grid methods. Nith the additional 
advantage of minimal memory requirements to store the coefficients of the 
governing equations as compared to the other methods, it is clearly the best 
for this problem. For the finest grid (162 x 162), the ratio of CPU times 
(best single-grid/best multi-grid) was 81 at Reynolds number of 100, but only 
14 at Reynolds number of 1000. Computer core memory requirements were 
1.50, 1.52, and 3.30 M Bytes for the multi-grid codes with point, line, and 
plane relaxations, respectively. 

Figure 5 shows some comparisons of the computed velocity profiles with 
those of Ghia et al . (ref. 12) at Reynolds numbers of 100 and 1000. The agree- 
ment is seen to be perfect, within the limits of accuracy of the plots. The 
predicted sizes and strengths of the primary, secondary, and tertiary vortices 
also compare quite well. Figure 6 shows the streamline plots at both Reynolds 
numbers . 

Table II shows the comparison of the CPU times and the number of itera- 
tions required for convergence for the backward-facing step flow problem. The 
comparisons are for 6 grid levels. The results show that, for single grids, 
the point-relaxation method exhibited the slowest convergence rate, and this 
became worse, the higher the grid aspect ratio. The coupled line- and plane- 
relaxation schemes had better convergence rates; albeit worse than the better 
decoupled methods. The deterioration of the convergence rate of the point- 
relaxation scheme with increased grid aspect ratio was even more pronounced in 
the multi-grid version. Typically, there was a four- to six-fold increase in 
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the CPU times for convergence when the grid aspect ratio was increased from 
3 to 6. Larger aspect ratios are required for higher Reynolds number computa- 
tions since the recirculation zone and the recovery length Increased with 
Reynolds number. Preliminary calculations showed that there was no Reynolds 
number dependence of convergence rates In this case. The multi-grid methods 
using line- and plane-relaxation schemes showed only modest increases in the 
required CPU times. The line-relaxation scheme appeared to have the best over- 
all performance. Brandt (ref. 11) did point out that point-relaxation schemes 
would not be efficient for equations with large differences in magnitudes of 
the directional coefficients, as would result from grids with high cell aspect 
ratio. The asymptotic convergence rates of point-relaxation schemes such as 
the Gauss-Seidel were shown to be inferior to those of line-relaxation schemes. 
The present results show that the use of the appropriate multi-grid method 
enabled the required CPU times, on the finest grid, to be reduced by factors 
of up to about 20 and 13 for L/H = 3 and 6, respectively, as compared with the 
best single-grid method. The factors are correspondingly lower for the coarser 
grids. In fact, significant advantages of the multi-grid methods on.y occurred 
for grids finer than 34 x 34. Below this, the best multi-grid method might 
converge slower than the best single-grid one. It should be noted, however, 
that the multi-grid concept always resulted in faster convergence, if the same 
relaxation scheme was utilized. Also, the linearity of the CPU times for con- 
vergence with the total number of grid nodes is again shown in the log-log 
plots in figures 7 and 8. A useful measure of the utility of the multi-grid 
technique is the ratio of the CPU times for best multi-grid method/ best 
single-grid method. Table III presents this ratio for the two test cases. 

Finally, the three multi-grid methods were applied to the lid-driven 
cavity flow at Reynolds number of 100, but at various aspect ratios (B/H). 

This test case combines aspects of the first two test cases; the flow reversal 
in the cavity with high cell aspect ratios. The computed results are similar 
to those of Pan and Acrivos (ref. 23), with a single central vortex at the low 
aspect ratios and multiple central vortices at higher ratios of 3 and 6. The 
convergence characteristics are presented in table IV. The convergence rates 
for the three calculation methods are similar for aspect ratio of 1 . On 
decreasing the aspect ratio to 1/6, through 1/3 all methods experienced much 
slower convergence rates. CPU times for convergence increased by factors of 
7, 34, and 18, for the calculation methods with plane-, line-, and point- 
relaxation schemes, respectively. Similarly, on increasing the aspect ratio 
to 6, CPU times for convergence increased by factors of 6 and 7 for the plane- 
and point-relaxation schemes, respectively. No convergence solution was 
obtained with the line-relaxation scheme. The residual norm decreased by two 
orders of magnitude and, then subsequently, oscillated about this value. Line 
relaxation is, of course, directionally sensitive. It is most effective along 
lines perpendicular to the main flow direction, sweeping from upstream to down- 
stream. In the present flow problems, the relaxation line would be mostly par- 
allel to the velocity vectors when the aspect ratio is 1/6, but perpendicular 

to them when the aspect ratio is 6. So the slowing down of the convergence 

rate in the former is explicable. But in the latter, the velocity vectors 
have reversed directions near the lid and on the opposite side of the cavity. 
Whereas, line relaxation is carried out from upstream to downstream with 
respect to the flow in the top half, the situation is reversed with respect to 
the flow in the bottom half. This latter situation is believed to be the 

cause of the oscillation and lack, of convergence in this case. The directional 

sensitivity of line-relaxation schemes make them unsuitable for use in general 
purpose codes. 
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Point- and plane-relaxation schemes do not have such directional sensitiv- 
ity, but the former Is more adversely affected at all aspect ratios differing 
from unity. Plane-relaxation schemes thus appear to be the safest compromise 
for a general-purpose multi-grid code. But there Is a penalty to pay for 
Increased computer core memory requirements. Point-relaxation schemes would, 
of course, be best If large cell aspect ratios can be avoided. 


CONCLUSION 

It has been shown that considerable improvements In convergence rates can 
be achieved In two-dimensional fluid flow computations on fine grids by 
employing the multi-grid concept. For nearly homogeneous difference equations, 
the multi-grid method with a point-relaxation scheme proposed by Vanka 
(ref. 13) gives the best performance in terms of both CPU times and computer 
memory requirements. However, this method Is inadequate for solving inhomoge- 
neous equations, resulting from large grid aspect ratios, or similar. In these 
cases, a line- or plane-relaxation scheme should be utilized with the multi- 
grid method, but the computer memory requirements would increase accordingly. 

In the particular cases presented here, plane relaxation appeared to be the 
most robust although it would not necessarily procure the fastest convergence 
1 n al 1 problems . 

The multi-grid methods showed a linear Increase In CPU times for conver- 
gence, as the computational grid was increased, whereas the single-grid 
methods exhibited a nearly quadratic Increase. Thus, for very fine-grid compu- 
tations multi-grid methods would always be better than any single-grid method. 
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TABLE I. - CPU TIMES AND NUMBER OF ITERATIONS FOR LAMINAR FLOW IN A LID-ORIVEN SQUARE CAVITY 

(a) Re = 100 


Gri d 

<NI x NJ) 

Si ngl e-gri d d 

Mul ti-gri d^ 

SIMPLE 

SIMPLEC 

PI SO 

SCGS 

CLSOR 

SIM 


CLSOR + FMG 

SIM + FMG 

7x7 
12 x 12 
22 x 22 
42 x 42 
82 x 82 
162 x 162 

0.4/29 
1 .4/55 
13.5/162 
184/548 
2559/1848 
N/A c 

0.5/42 
1.1/38 
5.5/61 
65/187 
885/625 
11 800/2000 

0.4/18 

1.1/26 

5.9/47 

70/42 

1010/472 

N/A c 

0.3/27 

1.9/91 

23.8/324 

254/853 

3240/2786 

N/A c 

0.1/23 

1.3/62 

16.4/194 

219/641 

2931/2128 

N/A c 

0.4/17 

2.0/54 

16/183 

255/617 

3600/1755 

N/A c 




0.8/19/35 

3.9/19/43 

13/14/34 

38/12/27 

145/11/23 

1.0/30/43 

4.1/20/44 

15/17/38 

55/15/34 

219/15/34 

2.8/16/69 

7.5/14/75 

16/15/38 

59/14/32 

245/11/41 


(b) Re = 1000 


7 x 7 

0.5/37 

0 6/53 

0.4/21 

0.4/34 

0.2/29 

0.5/22 



— 

12 x 12 

1.3/49 

1 .2/45 

1 .2/29 

1 . 5/68 

1.4/65 

1.5/43 

0.9/36/40 

1 .4/36/57 

4.6/25/110 

22 x 22 

8.3/97 

5.3/59 

4.7/36 

14.3/192 

14.3/168 

10/108 

4.6/20/49 

11.2/42/13 

16/32/157 

42 x 42 

81/244 

38/111 

34/69 

206/702 

145/420 

107/269 

26/23/68 

42/33/105 

44/28/129 

82 x 82 

950/687 

348/248 

373/182 

3045/2538 

1499/1006 

1492/728 

102/18/66 

131/22/82 

153/22/103 

162 x 162 

N/A c 

4496/750 

4957/549 

N/A c 

N/A c 

N/A c 

325/13/52 

315/12/47 

342/16/63 


d CPU time (sec) on S7881/number of iterations. 

b CPU time (sec) on S7881 /number of fine grid iterations/total work units. 
c Not available due to excessive CPU time. 


TABLE II. - CPU TIMES AND NUMBER OF ITERATIONS FOR LAMINAR FLOW OVER A BACKWARD FACING STEP 

(a) Re = 100, L/H = 3 


Gri d 

(NI x NJ) 

Single-gri d a 

Mul ti -grid b 

SIMPLE 

SIMPLEC 

PIS0 

SCGS 

CLSOR 

SIM 

SCGS + FMG 

CLSOR + FMG 

SIM 4 FMG 

6x6 
10 x 10 
18 x 18 
34 x 34 
66 x 66 
130 x 130 

0.4/35 

0.8/39 

5.0/88 

61/278 

812/922 

N/A c 

0.5/51 
0.9/46 
2 . 5/42 
26/120 
342/376 
4400/1180 

0.4/13 

0.7/19 

2.8/32 

32/99 

340/312 

N/A c 

0.1/36 

1.5/113 

15/305 

154/775 

N/A c 

N/A c 

0.2/15 

0.5/37 

6.5/122 

91/418 

900/1031 

N/A c 

0.1/11 

0.3/27 

4.5/93 

74/296 

760/1020 

N/A c 




0.8/31/55 

3.8/28/65 

16/22/68 

63/22/65 

214/21/56 

0.4/20/26 

2.0/23/24 

11/27/46 

55/30/55 

243/29/58 

0.4/28/34 

2.3/24/43 

12/23/56 

56/25/62 

257/26/67 


(b) Re = 500, L/H = 6 


6x6 
10 x 10 

0.4/35 

1.0/48 

0.5/52 
1 .0/51 

0.4/18 

0.8/24 

0.4/115 

4.2/314 

0.6/16 

0.7/30 

0.1/11 

0.4/33 

2.8/75/181 

0.4/24/29 

0.5/29/420 

18 x 18 

5.0/82 

3.6/59 

3.4/39 

50/1022 

4.0/68 

3.0/788 

17/78/274 

2.6/27/43 

3.1/33/657 

34 x 34 

74/338 

28/123 

40/123 

364/1846 

50/229 

94/354 

87/55/345 

18/37/725 

33/66/169 

66 x 66 

1000/1100 

371/404 

502/387 

N/A c 

900/1031 

964/1329 

358/51/365 

70/27/70 

95/22/110 

130 x 130 

N/A c 

4490/1200 

N/A c 

N/A c 

N/A c 

N/A c 

1200/41/292 

339/34/82 

318/25/90 


d CPU time (sec) on S7881/number of iterations. 

b CPU time (sec) on S7881/number of fine grid iterations/total work units. 
c Not available due to excessive CPU time. 
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TABLE III. - RATIOS OF CPU TIMES 
[Best single-grid/best multi-grid.] 



(3) SYMMETRIC GRID NODES. 



(b) CONTROL VOLUMES FOR THE I, Jth CELL. 


FIGURE 1. - STAGGERED GRID SYSTEM. 


TABLE IV. - PERFORMANCE OF MULTI-GRID METHODS IN LID- 
DRIVEN FLOW IN CAVITIES WITH DIVERSE ASPECT RATIOS 


Aspect 

Re = 100, 5-level multi-grid, 

66 x 66 grid 

(B/H) 

SCGS + FMG a 

CLSOR + FMG a 

SIM 4 FMG a 

1/6 

513/91/484 

953/150/1088 

237/122/276 

1/3 

129/35/122 

152/ 58/175 

95/ 47/115 

1 

29/11/37 

29/ 16/36 

32/ 13/27 

3 

56/27/54 

144/ 68/181 

64/ 34/65 

6 

189/97/179 

Osc i 1 1 ati ng , 

184/111/205 



no convergence 



a CPU times (sec) on S7881/number of fine grid itera- 
tions/total work units. 


U 0 



(3) LID-DRIVEN CAVITY FLOW. 



(b) FLOW OVER A BACKWARD -FACING STEP. 


FIGURE 2. - ILLUSTRATION OF TEST CASES. 
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CPU TIME, s CPU TIME 






<b) V-VELOCITY ALONG HORIZONTAL CENTRE LINE. 

FIGURE 5. - COMPARISON OF VELOCITY PROFILES: PRESENT 

CALCULATIONS. OOO CALCULATIONS OF GHIA ET AL. (12). 


18 






(b) Re = 100. 

FIGURE 6. - STREAMLINE PLOTS FOR L ID-DftIVEN CAVITY FLOWS SHOWING 
NORMALIZED STREAM FUNCTION (4>/U 0 h) CONTOURS. 
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